New Insights into the Volume Isotope Effect of Ice Ih from Polarizable Many-Body Potentials

The anomalous volume isotope effect (VIE) of ice Ih is calculated and analyzed based on the quasi-harmonic approximation to account for nuclear quantum effects in the Helmholtz free energy. While a lot of recently developed polarizable many-body potential functions give a normal VIE contrary to experimental results, we find that one of them, MB-pol, yields the anomalous VIE in good agreement with the most recent high-resolution neutron diffraction measurements—better than DFT calculations. The short-range three-body terms in the MB-pol function, which are fitted to CCSD(T) calculations, are found to have a surprisingly large influence. A vibrational mode group decomposition of the zero-point pressure together with a hitherto unconsidered benchmark value for the intramolecular stretching modes of H2O ice Ih obtained from Raman spectroscopy data unveils the reason for the VIE: a delicate competition between the latter and the librations.

N uclear quantum effects can manifest themselves quite prominently in macroscopic thermodynamic properties, like for example phase transition enthalpies, 1 negative thermal expansion (NTE), or density change at low temperatures upon substitution of a light by a heavier isotope. 2 The latter is called the volume isotope effect (VIE) and originates from the zeropoint energy of the lattice vibrations (phonons). Most materials show a normal VIE, which means substitution with heavier isotopes results in a smaller molar volume at temperatures approaching the absolute zero. A hand-waving rationalization in a classical picture is that the larger vibrational amplitude and concomitant volume ascribed to a lighter isotope compared to a heavier isotope, which both experience the same chemical interaction potential at the same temperature. In ice Ih, the most common form of solid water on earth, the VIE is anomalous, resulting in a smaller unit cell and thus molar volume of the H 2 O compared to the D 2 O isotopologue (about 0.1% up to 200 K 3,4 ). Despite its small magnitude, the effect has been very well quantified experimentally. Only recently, high-resolution neutron powder diffraction measurements 4 have reduced the uncertainties for the unit cell volumes of ice Ih compared to earlier work 3 over a wide temperature range and thus provide an excellent benchmark for atomistic interaction models that can be employed in computational studies.
Computational modeling of the VIE is very challenging. So far, water force fields ranging from simple fixed point charge up to sophisticated polarizable models have all predicted a normal VIE for ice Ih. 5−7 Density functional theory, on the other hand, is able to model the VIE of different ice phases. 7−9 However, even a qualitatively correct description depends very strongly on the computational settings, in particular the choice of the exchange-correlation functional. 7 Quantification of the VIE using embedded-fragment ab initio second-order many-body perturbation (MP2) theory has fared somewhat better. 10 But also here, the results are very sensitive to computational details like the basis set size and the embedding field. This also holds for the individual contributions of the different groups of phonon modes to the zero-point pressure, which are ultimately responsible for the VIE. 8,10 These contributions are commonly expressed in the form of mode-specific Gruneisen parameters and have not been benchmarked against experimental data. Consequently, a detailed understanding of the VIE's origin in terms of the competition of different contributions to the chemical interaction potential has been elusive so far.
In this work, we provide that understanding based on recently developed polarizable many-body potentials as interaction models. Building on our recent studies, 11,12 we employ the quasi-harmonic approximation to account for nuclear quantum effects in the Helmholtz free energy by means of extensive and high-precision phonon calculations. We find that the MB-pol interaction model, whose short-range part is rooted in coupled-cluster calculations, yields the anomalous VIE of ice Ih in better agreement with the experimental reference value than DFT calculations with the PBE functional.
A decomposition of the zero-point pressure into contributions from different vibrational mode groups together with a hitherto unconsidered benchmark value, which we obtain from Raman spectroscopy, 13 allows us to scrutinize this further. According to the MB-pol total energy partitioning, the delicate competition between the librational and intramolecular stretching modes driven by a surprisingly large influence of short-range three-body effects is responsible for the anomalous VIE of ice Ih.
For the quantification of the VIE, this work employs the quasi-harmonic approximation (QHA), which has been used successfully for the same purpose in the past. [7][8][9][10]14 According to the QHA, the Helmholtz free energy of the ice crystal is given by and is conveniently evaluated per molecule. U is the internal energy that describes the interaction between molecules in the crystal. ω i are the vibrational modes which determine the second (zero-point energy E ZP ) and temperature-dependent third term of F. For the sake of simplicity, they are denoted by a collective index i that stands for both wave vector and band indices of the corresponding phonon modes. The internal energy U and the vibrational modes ω i depend on the unit cell volume V, so that the minimum of F with respect to V at a given temperature is generally different for H 2  In this work, calculations with different interaction models are performed, employing and extending the Atomic Simulation Environment (ASE) 15 for interfacing their respective implementations. This includes the fixed-pointcharge-based force fields q-TIP4P/F 16,17 (as available in LAMMPS 18 ), the polarizable force fields AMOEBA14 19,20 (as implemented in TINKER 21 ), SCME/f, 12,22 and MB-pol (as implemented in the MBX package). 23−25 All-electron density functional theory (DFT) calculations with the PBE exchangecorrelation functional 26 are performed with the FHI-aims code, 27,28 using the same high-accuracy settings thoroughly verified 29 and employed 11,30 for ice Ih in previous work (see the Supporting Information for details). The DFT calculations mimic proton disorder with a simulation cell containing 12 molecules. 31 A simulation box containing 96 molecules, or some larger supercell, have been used for the force field interaction models to ensure the same level of convergence for V 0 (±0.01 Å 3 per molecule 29 ). For all interaction models, V 0 is calculated with ASE as in our earlier work by a combined optimization of the cell vectors and the molecular degrees of freedom preserving the space group of the lattice (as defined by the oxygen atoms) 11 with a maximum force threshold of 1.0 × 10 −3 eV Å −1 .
A continuous representation of U(V) is obtained by leastsquares fitting to the Rose−Vinet 32 equation of state. Isotropic contraction and expansion of V 0 by ±4% yields 11 structures for each interaction model, for which again all molecular degrees of freedoms have been relaxed. Phonon calculations have been performed for all of these structures with the PHONOPY code, 33 using a finite displacement 34 of 0.02 Å in 3 × 3 × 3 supercells of the original simulation cell. The Brillouin zone has been sampled by 30 × 30 × 30 and 10 × 10 × 10 grids of phonon wave vectors in the 12 molecule and 96 molecule simulation cells, respectively. The implementation of the QHA in PHONOPY then yields a continuous representation of the volume-dependent second and third terms in eq 1 and thus , 35 and the phonon mode-dependent Gruneisen parameters γ i . Convergence checks for the VIE can be found in the Supporting Information. Table 1 compiles experimental data and results from calculations at both 0 and 10 K to confirm that this has no effect on any of the numbers presented in Table 1. As demonstrated only recently, errors related to the treatment of core and valence electrons in different DFT codes can be sizable for the calculation of energy−volume curves U(V) 36 and could thus significantly affect results for the VIE. This is the most likely reason for the difference of 0.07% between earlier DFT results obtained with the same exchange-correlation functional (PBE). 7,14 Our own PBE calculations eliminate this source of error and perfectly The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter reproduce the value for V 0 obtained in earlier work. 29,30 In combination with our meticulous convergence tests for the phonon calculations with respect to the VIE (see the Supporting Information), we can therefore confirm without any further doubts the conclusions from earlier work, 7,14 namely that the PBE functional reproduces the experimentally observed anomalous VIE but overestimates it. Likewise, our computational setup also confirms that the q-TIP4P/F force field, which is based on fixed point charges, yields a normal VIE for ice Ih. 7 The same holds for the two polarizable force fields AMOEBA14 and the recently established SCME/f. The results for V H O 2 and also for V D O 2 obtained with SCME/f show the best agreement with the experimental values. SCME/f is followed by MB-pol, which also yields and consequently a correct description of the anomalous VIE. The absolute value of 0.14% is even in much better agreement with the experimental data than (our) PBE results. This remarkable result also provides the opportunity to better understand what contributions to the chemical bonding in the ice Ih crystal are responsible for the VIE. Among all the force fields considered here, MB-pol is the only one that explicitly accounts for shortrange interactions involving triples of water molecules, which have been parametrized to quantum-chemical CCSD(T) calculations. 23 Indeed, omitting these terms (MB-pol w/o 3B in Table 1) yields a normal VIE.
At first glance, the strong influence of these terms is surprising because they do not constitute a large contribution to the cohesive (lattice) energy according to the decomposition of the internal energy at the equilibrium volume U(V 0 ) in the MB-pol force field. 38 This does not change when moving away from V 0 as shown in Figure 1. Compression of the ice Ih lattice leads to an increase of the repulsive shortrange interactions between pairs of water molecules (two-body terms), which is almost compensated by the increase of the long-range electrostatic and dispersion contributions. The opposite holds when expanding the volume. The one-body (i.e., deformation of individual water molecules) and shortrange three-body terms play hardly any role.  Figures 2a and 2b, respectively. The experimental data for H 2 O and D 2 O ice Ih shows a negative slope for T ≤ 60 K. This negative thermal expansion (NTE) has been modeled successfully before 10,39 and is also reproduced by all methods considered here, except for SCME/f. This leads to a small offset in the relative volume change for H 2 O between SCME/f and the experimental data at higher temperatures, which remains almost constant. Apart from that, SCME/f captures the shape of the experimental curve for  (2) Positive (negative) values for P ZP yield expansion (contraction) of the volume due to zero-point energy effects. Figures 3a  and 3b show that all methods yield a positive total zero-point pressure for both H 2 O or D 2 O, respectively, as to be expected according to the results compiled in Table 1. For an anomalous (normal) VIE (at 0 K), the total P ZP needs to be smaller (larger) for H 2 O than for D 2 O. The corresponding differences are shown in Figure 3c, and indeed only PBE and MB-pol yield . Figure 3 also shows a decomposition of the zero-point pressure into contributions from the five different vibrational mode groups characterized by hydrogen-bond bending (HB) and stretching (HS), librations (L), intramolecular bending (B), and stretching (S). Unlike the differences of the total P ZP , which unfortunately cannot be measured directly, the differences between the contributions from the mode groups vary much more when considering different interaction models. The P ZP contribution from the HB and HS groups is hardly affected by H−D substitution (see Figure 1. Lattice energy E lat and its contributions according to the total energy decomposition of the MB-pol force field 23,24,37 (all in eV per molecule). The intramolecular (1B) as well as intermolecular short-range two-body (2B), short-range three-body (3B), long-range electrostatic (elec), and dispersion (disp) contributions add up to E lat (total). Violet, red, and blue bars depict equilibrium (V 0 ) and isotropically compressed (0.96V 0 ) and expanded (1.04V 0 ) lattice configurations, respectively, as encountered during the phonon calculations for the determination of the VIE according to the QHA.
The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter Figure 3c), which is not surprising because both of these mode groups involve the frustrated translation of entire H 2 O and D 2 O molecules. Likewise, all methods suggest that the B modes contribute very little to the VIE and that it is the delicate balance between the expansive P ZP of the L modes (frustrated rotations) and the contractive P ZP of the S modes, which predominantly determine the sign of the zero-point pressures difference between both isotopologues. Salim et al. 10 have already pointed that the subtle interplay of the P ZP contributions from different mode groups makes it additionally challenging to determine whether a particular interaction model captures the VIE correctly for the right reason. We have already noted in our earlier work 11 that the measurements by Minceva-Sukarova et al. 13 of the pressure dependence of the Raman peak for the S mode group in H 2 O ice Ih at 246 K play a key role in this context. As further detailed in the Supporting Information, this allows us to obtain a good estimate for the contribution by the S mode group P ZP S −0.548(51) GPa, which is based on experimental data alone. Figure 3a includes this value as a black horizontal line. Considering that error estimates are lower bounds, MB-pol almost reproduces this value exactly (−0.66 GPa) and clearly comes much closer than PBE (−1.32 GPa) and MB-pol w/o 3B (−0.33 GPa). This confirms that short-range three-body effects play indeed a very important role for the correct atomistic description of the VIE.
Unfortunately, Minceva-Sukarova et al. 13 have not measured the S mode group frequency shift for the D 2 O isotopologue of ice Ih. Because all interaction models suggest a strongly localized character of the S modes, the outcome of such a measurement can be estimated based on the reduced masses associated with the O−H and O−D bonds as P 0.728 ZP S −0.399 GPa (see the Supporting Information), which is consequently again in excellent agreement with the MB-pol value in Figure 3b (−0.47 GPa). While it would be good to see this value confirmed in future experiments, the outstanding performance of MB-pol and the concomitant understanding would be better scrutinized by experimental data for the L modes. Ideally, such data could also be measured at temperatures much lower than 246 K to avoid systematic errors related to mode softening effects in future comparisons of experiments and theory.
In summary, this study provides new insights into the volume isotope effect of ice Ih based on calculations within the quasi-harmonic approximation by employing a variety of different interaction models. Numerically precise all-electron = Table 1 for the same interaction models except AMOEBA14. (c) Resulting temperature dependence of the volume isotope effect